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Massive bigravity, a theoretically consistent modification of general relativity with an additional 
dynamical rank two tensor, successfully describes the observed accelerated expansion of the Uni¬ 
verse without a cosmological constant. Recent analyses of perturbations around a cosmological 
background have revealed power law instabilities in both scalar and tensor perturbations, moti¬ 
vating an analysis of the initial conditions, evolution, and cosmological observables to determine 
the viability of these theories. In this paper we focus on the tensor sector, and study a primor¬ 
dial stochastic gravitational wave background in massive bigravity. The phenomenology can differ 
from standard General Relativity due to non-trivial mixing between the two linearized tensor fluc¬ 
tuations in the theory, only one of which couples to matter. We study perturbations about two 
classes of cosmological solutions in bigravity, computing the tensor contribution to the temperature 
anisotropies in the Cosmic Microwave Background radiation and the present stochastic gravitational 
wave background. The result is strongly dependent on the choice of cosmological background and 
initial conditions. One class of background solution generically displaying tremendous growth in the 
amplitude of large-wavelength gravitational waves, while the other remains observationally indistin¬ 
guishable from standard General Relativity for a wide variety of initial conditions. We analyze the 
initial conditions for tensor modes expected in an inflationary cosmology, finding again that there is 
a strong dependence on the assumed background. For one choice of background, the semi-classical 
theory is beyond the perturbative regime. For the other choice, inflation generically yields initial 
conditions that, when evolved, give rise to a stochastic background observationally indistinguishable 
from standard General Relativity. 
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I. INTRODUCTION 


Extending general relativity (GR) by giving a mass to the graviton has proven to be a difficult endeavor, dating 
back to first attempts by Fierz and Pauli in 1939 [lj. As a result of recent work in constructing a fully non-linear 
ghost-free theory, dRGT massive gravity [2], interest in this field has been on the rise (see [3] for a review). An 
important yet curious feature of this theory is the need to introduce a fixed non-dynamical tensor field, in 

addition to the metric g^ describing spacetime. Promoting this tensor to a dynamical field, a ghost-free theory of 
massive bigravity [4j has recently been introduced. With two dynamical metrics, the choice of matter coupling is a 
non-trivial issue [5j; we will consider couplings of matter to one metric only. In this scenario, one metric g^ describes 
our spacetime, and another one, f^ v , is part of a “dark” gravitational sector. 

Much of the interest in massive gravity is motivated by the puzzle of cosmic acceleration and dark energy. In the 
standard cosmological model, this accelerated expansion is assumed to be due to a cosmological constant. However, 
the extreme fine tuning of the cosmological constant has led physicists to pursue alternative explanations for cosmic 
acceleration. Viable homogenous and isotropic cosmological solutions exist in bigravity which can describe our uni¬ 
verse, including acceleration, without a cosmological constant [ME]- These are usually referred to as self-accelerating 
background solutions. Beyond the background level, investigations have recently been underway to analyze pertur¬ 
bations in bigravity[T2HT7]- These analyses reveal that there is a particular class of stable solutions, but all others 
are plagued by an exponential instability in the scalar sector in the early universe j!2] . 


The transverse traceless fluctuations of each of the two dynamical metrics in bigravity interact, altering the prop¬ 
agation of gravitational waves as compared to GR. The most general set of linearized equations of motion for the 
visible sector tensor modes h g and the dark sector tensor modes hf is given by 


where 


D 2 ■ h + m 2 (x, r) • h = 0 




m 2 (x, r) 


m|(x,r) mj (x,r)\ 
m s/( x ’ r ) TO s( x ’ r ) ) ’ 


(1) 

( 2 ) 


where V 2 is the covariant derivative defined with respect to / M „ and V 2 is the covariant derivative defined with re¬ 
spect to g M „. There is no off-diagonal term in the differential operator D 2 due to the absence of consistent derivative 
couplings between the two metrics. The mass matrix m 2 is in general not diagonal or symmetric. The differential 
operator D 2 and mass matrix m 2 are generally not simultaneously diagonalizable, leading to mixing between the 
visible and dark sector tensors. Because the visible and dark sector tensor fluctuations evolve in a different (spacetime 
dependent) background and possess a different (spacetime dependent) mass, their propagation speeds can in general 
be different and spacetime dependent. The modifications in the propagation speed of gravitational waves and mixing 
between the visible and dark sector tensors generically present in bigravity can have important implications for both 
astrophysical [TM2T] and primordial [221427] gravitational waves. 


In this paper, we consider a stochastic background of primordial gravitational waves in massive bigravity, studying 
the impact on primordial and present day gravitational wave observables by computing the tensor contribution to 
the temperature anisotropies of the Cosmic Microwave Background radiation (CMB) and the power spectrum of 
the present day stochastic gravitational wave background. We examine whether these observables can be used to 
distinguish between bigravity and GR, serving as a possible test for gravity on cosmological scales. We consider two 
classes of cosmological background solutions for the dark sector metric f gv . For one class, the tensor perturbations 
match those of GR very closely. However, for the second class of background solution, a power-law instability in the 
tensor sector leads to a strong growth in the amplitude of gravitational waves at late times. In the absence of special 
initial conditions, these solutions are in conflict with observations. We argue that these conditions may be naturally 
produced in the standard inflationary picture. In addition, we explore to what extent this fine tuning is present in 
solutions that are not self accelerating. 


During the preparation of this manuscript, the analysis of Cusin et. al. nz! appeared; additionally the analysis 
of Amendola et. al. [27l appeared. There is significant overlap in our discussion of the evolution of perturbations. 
However, we make substantially different assumptions about the initial conditions, which we motivate from a detailed 
discussion of bigravity in the context of inflationary cosmology. This analysis addresses an outstanding question 
regarding the viability of massive bigravity first outlined in Cusin et. al. We explicitly comment on the relevant 
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differences in our assumptions throughout the manuscript. 

The paper is organized as follows: In section [II] we give the bigravity background equations and the two classes of 
solutions under consideration. Section |III| analyzes the gravitational waves in each branch. In the following sections 
we compute gravitational wave observables: section |IV| gives the CMB Tensor Power Spectrum in bigravity, and 
section |V| gives the present day stochastic gravitational wave background. Section pVT] explores the initial conditions 
as predicted by inflation, followed by a discussion in | VII[ 


II. BIGRAVITY COSMOLOGY 


We will consider the massive bigravity model proposed by Hassan and Rosen |¥j: 


S = -~Y [ d A x^gR(g) - ~y f d A x^fR{f)+m 2 M 2 f Priori ( \J~ 9 */) + ^matter (3) 

■' ■’ n =0 

where the two dynamical metrics are and / M „ and their respective Planck masses are M g and My, while m is 
the mass scale associated with the graviton mass matrix. The interaction term between g and / contains a linear 
combination of the symmetric polynomials e n (X) which are defined by 


eoPO =1 

(4) 

ei(X) =[X) 

(5) 

e 2 (X) =i([X] 2 - [X 2 ]) 

(6) 

e 3 (X)=i([X] 3 -3[X][X 2 ]+2[X 3 ]) 

(7) 

e±{X) =det(X) 

(8) 


given in terms of the tensor X£ = y g^fav where square brackets denote the trace. The dimensionless coefficients 

Pn are free parameters of the theory. Further, we will consider a singly-coupled theory in which Scatter contains only 
couplings of matter to g^. In this case, g„„ is considered the standard physical metric while f^ is a new dynamical 
tensor field. This singly-coupled theory was shown to be free of the Boulware-Deser ghost [5]. The resulting quantum 
corrections (at one loop) are nothing other than the standard cosmological constant which does not detune the special 
structure of the graviton potential and is thus harmless [5j. By rescaling the dark metric and the free parameters as 
follows, 

= (Mf/M g )% v PI = (M g /Mf) n p n , (9) 

we can rewrite the action so that the redundant scale Mf is absent: 


S = Mg 


j d 4 Xy/=gR{g) j d A xy / ^fR(f) + m 2 J + 




( 10 ) 


It has been shown that the cosmological expansion and spherically-symmetric solutions to this theory give viable 
alternatives to general relativity (ME! El] , at the background level. Let us briefly examine the background equations 
of motion and solutions. 


We start by making an FRW ansatz for the metrics g^ and 

ds 2 g =a 2 (r)(— dr 2 + dxidx 1 ) (11) 

ds 2 =b 2 (t)[—c 2 (r)dT 2 + dxidx l ] (12) 

where r represents conformal time, a and b are the scale factors corresponding to the g and / metric respectively, and 
c is the lapse function for the / metric. Note that the assumption that is FRW is not the most general choice for 
the metric, and it could have more complicated dynamics. 
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The Bianchi identities imply the following relation: 


H f ba r 

c = — J- = — = 1 -t-—. 

ab 


n 


rH 


(13) 


where 7~L = a/a is the conformal Hubble function for the g metric, Hf = 6/6 is the conformal Hubble function for the 
/ metric, and r = 6/a the ratio of scale factors. For simplicity we perform the following rescaling 


b =?ZL/r 

yn — tt 2 yn 


P = 


UgHy 


v- n 

H -w 0 


(14) 


where p = p m + p r is the dimensionless energy density of all matter and radiation components, and we measure all 
times and lengths in terms of Ho. With these definitions in hand, variation of the action with respect to g^, and 
inserting (11) and (12), gives the Friedmann equations 


3'H 2 = a z {p + p mg ) 
a 3 dp 


(15) 

2 U + 'H 2 = a 2 p + -r-j- + a? ((3o + 0ir(2 + c) + 0 2 r 2 (l + 2c) + /3 3 r 3 c )) (16) 

3 da v y 

where p mg = (/3 0 + 3/3i r + 3/3 2 r 2 + ^ 3 r 3 ) (17) 

where p mg is an effective massive-gravity energy density. We also have the background equations for the / metric: 

ji 


3 H 2 = — (/3i + 3/3 2 r + 3/3 3 r 2 + /d 4 r 3 ) 

a 2 

Tti + H 2 c = — (0i + 0 2 r{2 + c) + /3 3 r 2 (l + 2c) + 0 4 r 3 c) 

The energy densities follow the usual conservation laws 

p m + 3Hp m = 0 and p r + +Hp r = 0 , 

giving rise to their solutions in terms of a and the present day density parameters = p®/(3 HqM 2 ) = p//3: 


Pm — 


o° 

^ L m 

3a 3 


and p r = 


3 a 4 


(18) 

(19) 

( 20 ) 

( 21 ) 


The above equations can be combined to form convenient equations for the dynamic variables r and a: 


0 = pm + Pr ~ ~ (0i + 3/3 2 r + 30 3 i ’ 2 + 0 4 r 3 ) + 0 O + 30ir + 30 2 r 2 + 0 3 r 3 

1 dr _ 3 0 3 r 4 + (3/3 2 - 0i)r 3 + 3(0i - 0 3 )r 2 + (0 O - 3 0 2 )r ~ 0i - H/q~ 4 r 

r da a 30 3 r 4 + 2(3/3 2 - 0 4 )r 3 + 3(0i - 0 3 )r 2 + 0i 


The procedure now is as follows: 

1. Fix the free parameters: fi°, 0 O , 0i, 0 2 , 03, 


2. Use (22) evaluated today to find r 0 , the present value of r. 


3. Starting with initial condition ro, evolve (23) to find r(a) for all time. 


4. Using r(a), solve for the Hubble parameter H(a) using (18). 


( 22 ) 

(23) 


5. It is now possible to find a(r) using H = a/a, and thus can also find r(r) and H(t) and c(r) 

In this work we will focus on two distinct types of background solutions, following the notation in [16] . A plot of 
r(a) in each branch is shown in FIG. [TJ 
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Expanding Branch Bouncing Branch 




FIG. 1: The evolution of the ratio of scale factors r = bja for the both branches described in the text. 


1. The Expanding Branch: In this branch of solution, both metrics g and / expand with time. This is also 
known as the finite branch since the ratio r = b/a evolves from zero at early times to a finite value. Within 
this branch, there is a proposed minimal model HE HU in which only Al A 0- At the background level, it was 
shown that this theory can be compatible with expansion histories, but remains distinct from GR with testably 
different observables HE HU. The issue with this branch is an exponential perturbative instability in the scalar 
sector, previously noted in the literature pi muni hu. We therefore focus on this branch as an example of how 
phenomenologically different the results for tensors can be under different assumptions about the background 
solution. We will work within the minimal model, fixing p 1 = 1.38, which was found from a best fit analysis at 
the background level in |B]. In this case, an analytic solution exists for r(a) and is given by 


r(a) 


-3 an° m - 3fi° + ^3(4a 8 /3? + 3afi°, 2 + (kA2°/2° + 3fi0 2 ) 

6a 4 /3i 


(24) 


2. The Bouncing Branch: This is a more exotic option in which the physical metric g expands in time while the 
dark sector metric / experiences a bounce. At the bounce point, /oo = 0, thus, /A, 1 diverges, but there is no 
divergence in the physical sector [Mini HU- In addition, well-defined and stable solutions for the background 
and linear perturbations exist through this point jl6| indicating that this divergence is likely of a mathematical 
rather than physical nature. The perturbative instability of the expanding branch is absent in this sector. This 
branch is also called the infinite branch since r evolves from infinity at early times to a finite value at late times. 
Here we are required to set A) = Pi = P 3 = 0, and for the remaining parameters, we use the best-fit values 
Pi = 0.48 and At = 0.94 found by fitting growth histories and type la supernovae Wi¬ 


lli. TENSOR PERTURBATIONS 


Here we give the equations of motion for the transverse traceless tensor modes /iJA and hjjj corresponding to 
the metrics g Ml/ and f respectively. To compute observables, we will be most interested in the perturbations 
corresponding to the physical metric, hJJj, since these are the ones that couple to matter. The tensor equations of 
motion were first analyzed in |13j . A thorough analysis of scalar, vector, and tensor perturbations was performed in 
m- In addition, further analysis of tensor perturbations appeared in Refs. m and m- The tensor perturbation 
equations of motion in momentum space are 

hg + 2'Hhg + k 2 h g + a 2 \(h g -hf) = 0 (25) 

hf + c 2 k 2 hf — ^ C (fe g — hf) = 0 (26) 


2 [%+-)- 


where superscripts and subscript indices have been dropped for simplicity. In addition, the time-dependent function 
A is defined as 


A = P 3 cr 3 + P 2 (c + l)r 2 + Au 


(27) 
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which simplifies to A = f3\r in either branch under consideration. These equations are satisfied separately for each 
polarization; the polarizations do not mix. 


Initial conditions for h g and hf at some initial time 77 are required to obtain solutions. In the absence of a theory of 
initial conditions, we should consider general initial data ht g j){ri) and /i( g j)(rj). However, for standard inflationary 
cosmology in GR, tensor modes free ze in once their physical wavelength becomes comparable to the primordial horizon 


size, motivating h g (ji) = 0. In Sec. VI we compute the initial conditions expected for inflationary cosmology in the 


context of bigravity, finding that ^( s ./) ( 77 ) = 0 is an appropriate choice. Note that this assumption differs from Cusin 

et. ah, who consider initial data with hf(Ti) 7 ^ 0. This was motivated by the presence of a growing mode for hf , 
which if excited, would dominate the evolution. Our choice of initial data initially sets this growing mode to zero, 
and as shown below, this leads to a different growth history for hf at late times. 


1. Expanding Branch 


In the expanding branch, the factor [2 (H + 7 ) — j] is always positive, and causes significant damping of the dark 
sector tensor perturbation hf. In addition, the factor a 2 A is small in this branch, a 2 A < 0.3 for all time. Therefore, 
unless the initial amplitude of hf is very large compared with h g , the mixing term a 2 \{h g — hf) does not significantly 
alter the behaviour of the physical tensor perturbation h g . For equal initial amplitude, our numerical solution for 
the tensor modes in this branch match closely with those of pure GR. This was confirmed for modes ranging from 
k = O.lHo to k = 10 5 Hq. Refer to FIG.kyfor the numerical solutions of this branch as compared to the standard GR 
gravitational waves with 77 = 10 — .TH 


If the relative amplitude of the dark sector tensor mode was decreased, hf{ri) < h g {ri ), this would only drive 
the physical tensor modes closer to those of GR. But if hf had a much higher initial amplitude relative to h g , the 
mixing term a 2 Xhf in (26) can become dominant for some time. However, even if we set 11 /( 77 ) h g (ji ), hf decays 
so dramatically that there is very little effect on h g . The bottom right plot in FIG. [2] shows how the dark sector 
perturbation decays very quickly, even when starting with a much larger amplitude. 


2. Bouncing Branch 


In the bouncing branch, the oscillations of hf are anti-damped for r < 77 , and damped for r > 77 , where 77 , is the 
bounce time corresponding to when c(r) = 0. This happens at relatively late times, around z\> ~ 0.6. Until the 
bounce occurs, while c < 0, hf experiences enormous growth, then starts to decay after 77 ,. Through the mixing term 
in (25), this growth in hf translates into growth in the physical mode h g , leading to oscillations that grow at late 
times. In FIG. [3j we show the evolution of the physical and dark sector tensors for equal amplitude initial conditions 
with Ti = 10 ~ 6 Hq 1 ; the large deviation from GR at late times is clear. 


The growth in h g (r) at late times is fc-dependent. Empirically, we find a falloff proportional to 1/fc 2 . This 
amplification is also dependent on the initial time 77 , which physically relates to the reheating temperature T^eat = T). 
This can be understood by examining the solution for h / within the radiation dominated era. On super-horizon scales 
there is an exact solution m given by 

hf=Ci+c 2 T 3 (28) 


so there is both a constant mode and a growing mode. Our choice of initial condition hf = 0 selects the constant 
mode (ie. c 2 = 0), so naively one might think that h f sh ould not grow at all in the radiation dominated era, regardless 
of the initial time 77 . However, the above solution (28) is only an approximate solution on super-horizon scales, not 
valid for k 7 ^ 0. For k / 0 we expect to depart from the the constant mode solution on a timescale of r ~ 1/cfc, at 
which point the growing mode will completely dominate. The earlier the initial time is set, the more time hf has to 
grow, ultimately driving more growth in h g . We find that the growth is inversely proportional to the initial time: 


1 Although this is choice for Ti, corresponding to reheat temperature T, = T eq [a eq /a(r,)] = To(ao/a e q) 2 [a e q/a(ri)] = 0.07 GeV, is not 
entirely plausible, it is still deep within the radiation era and is a practical choice for our numerical analysis. We will be able to 
extrapolate to earlier Ti when necessary using scaling properties of the solutions defined in (|29[). 
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Expanding Branch k = 100H 0 Expanding Branch k=100H 0 




Expanding Branch k = 100 H 0 



Expanding Branch k = 100 H 0 



FIG. 2: Top left: the solution for h g (r) at k = IOOHq for the expanding branch (green) versus GR (blue). The two solutions are 
essentially indistinguishable. Top right: the solution for hf(r ) at k = lOOHo for the expanding branch. Bottom left: The difference 
between the expanding branch solution and the standard GR solution, shown to be less than 5 X 10 -5 for k = 100i/o- Bottom right: 
the solution for hf(r) at k = lOOHo for the expanding branch, with hf(ri)/h g (ri) = 10 6 . Even in this case of hf{ji ) h g (ri ), the decay 
of hf is so fast that it does not cause any alteration to the physical tensor perturbations. 


h g (ro) oc 1 /i~i (xTi, valid for all initial times in the radiation dominated era. 

In addition to this, the growth is also proportional to the initial condition for hf. However, for a small enough 
value of hf(ri), the solution for h g {r) does not scale. For our choice of initial time r, = lO -6 #,^ 1 (corresponding to a 
reheat temperature of T* = 0.07 GeV) we find that for /i/(rj)//i s (r,) < 10 —9 , the solution for h g is indistinguishable 
from its solution with hf{ri) = 0 which agrees very closely with the pure GR solution. To extrapolate this result to 
a more reasonable reheat temperature, say Tj = 10 10 GeV, we must consider that hf grows proportionally to T t . We 
conclude that for Tj = 10 10 GeV, we need a suppression of hf{ri)/h g {ri) < 10~ 20 to obtain solutions that agree with 
those of GR. At this threshold, scaling down hffa) further will not result in any significant change. It is evident 
that to control the large growth in the bouncing branch, seeking gravitational waves that do not substantially deviate 
from those of GR, we require detailed knowledge about the mechanism by which they were produced. This will be 
explored further in Sec. [VT| 

From a purely phenomenological standpoint, one might also be interested in varying the graviton mass. This reduces 
the strength of the mixing term, which has the form oc m 2 /3*a 2 r(/i ff — hf)/H q (reinserting factors of m and Hq from 
©)• However, lowering m means weakening the influence of the bigravity interaction term in |3]) to the point at 
which it can no longer yield a viable background cosmology, so a cosmological constant must be reintroduced. This 
undercuts the strongest theoretical motivation for bigravity, but mild variations could be of interest in constraining 
the parameters of the theory. We find that h g (ro)/h g {ri) oc m. Just like for hffc), there is also a critical value of m 
for which this scaling relation no longer holds. Putting this all together, we find the following general relationship for 
the growth: 


h g (r 0 ) mhf(Tj) 
hg{Ti) ^ Tjfc 2 


mh f (Ti)H(Ti) mhf(Ti)Ti 
oc-oc 


k 2 


k 2 


M 7 *) 

hgi^y 


for 


m/H 0 > A clit 


(29) 



























Bouncing Branch 


Bouncing Branch 


k = 100 H 0 



k = 100 H 0 



FIG. 3: Top left: the solution for h g (r) at k = 100Ho for the bouncing branch (red) versus GR (blue). Top right: the solution for 
h f (r) at k = 100 Ho for the bouncing branch. 


k-dependent growth in bouncing branch 



r,-dependent growth in bouncing branch 



10" 4 0.001 0.010 0.100 
T H 0 


h/[r,(-dependent growth in bouncing branch 



m-dependent growth in bouncing branch 

100 ----—-- 

1 

E 001 

o> 
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10' 4 

icr 6 

10 -4 0.001 0.010 0.100 1 10 
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FIG. 4: The growth in the physical sector is dependent on k , Ti, hf(ri) and m < |29| ). We vary one parameter individually per plot to show 
how the solutions change, fixing all but one of = 10 —6 iif c j _1 (T* ~ 0.07 GeV), k = 100Ho, m = Ho, and h( g ,f)(ji) = 1. Top left: The 
gravitational waves in the bouncing branch grow like 1/k 2 . Here we plot k = 10 p ' 2 Ho with p ranging from 2 to 10 from top to bottom. 
Top right: The tensor modes grow like 1 /n. Here we plot Ti = 10 —8+p / 2 H ( ^~ 1 with p ranging from 0 to 8 from top to bottom. Bottom 
left: The gravitational waves grow proportional to hf(ri). Here we plot hf(ri) = 10~ p ' 2 with p ranging from 2 to 11 from top to 
bottom. Bottom right: The growth in the bouncing branch grow proportional to m (requires re-inserting A into the theory). Here we 
plot m = 10—P/ 2 Ho with p ranging from 2 to 11 from top to bottom. 
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where A c r it is the critical value at which the solution no longer scales (A cr it = 10 _2 ° for Tj = 10 10 GeV). Note that 
this scaling is different than that found in Cusin et. ah, which is due to the difference in choices for initial data. These 
growth dependencies are displayed in FIG. [4] which shows the solutions for various k, hffa), and to. Notice that 
the variations in growth of h g are similar for hf{ri) and ?n, which comes from the fact that both of these effects act 
to alter the coupling term in (251. Although the effect on growth in the physical sector appears the same, altering m 
introduces no change in growth in the dark sector. 


As shown in Ref. m , the generalized Higuchi bound in the tensor sector on an FRW background is satisfied for the 
expanding branch, but violated for the bouncing branch. This is essentially due to the behaviour of the lapse function, 
c, which takes on negative values in the bouncing branch, causing the kinetic term for the tensor perturbations to 
become negative. A violation of the Higuichi bound in an FRW background leads to a power law instability that is 
evident in our growing numerical solutions. This instability is not necessarily detrimental to the theory so long as 
these unstable modes are sufficiently suppressed, an important point which we investigate in section [VI] 


IV. CMB TENSOR POWER SPECTRUM 


In this section we use the solutions of (251 and (26) to compute the predicted bigravity CMB temperature tensor 
power spectrum (see m and [33] for recent computations of CMB spectra in bigravity). We begin by using the 
physical tensor mode h g (k , r) to compute the Ith photon moment due to tensor perturbations (assuming instantaneous 
recombination) using 


©F = J dr h g {k,T)ji[k(T 0 - r)]] (30) 

where r 0 is the conformal time today, r* is the conformal time at the time of last scattering, and ji is the spherical 
Bessel function. The tensor contribution to the temperature anisotropies is entirely due to the Integrated Sachs Wolfe 
effect. The evolution of the visible sector tensors can be substantially different than GR, leading to potentially visible 
differences in the spectrum of temperature anisotropies. For example, in the bouncing branch, from FIG. [4] it can be 
seen that there is additional time dependence on super horizon scales, and late-time growth of tensor modes; both of 
these effects will alter the temperature anisotropies. 


Once the photon moments are found, the angular power spectrum can be found from: 

ef 


,rp (l — 1)(Z + l)(l + 2) 


C{ = 


1 

dk - 
k 


«T 


+ 2 


U (l+ 2 ) 


(21 — 1) (2Z + 1) (21 — 1)(2Z + 3) (21 + l)(2l + 3) 


(31) 


We solve for the photon moments and angular power spectra numerically. 

For the expanding branch, the Cff’s are approximately the same as in GR, which is expected given the agreement 
of h g (r) with /igr(t)- However, in the bouncing branch, there is a drastic difference in the power spectrum. 


Starting with the measured value of the scalar quadrupole Cj GR ~ 1000/rA' 2 


and saturating the bound on the 
20CfoAT 2 . In d30b the 


tensor-to-scalar ratio of r = Cj/Cf < 0.2, the expected value of tensor quadrupole is Cj GR 
integral runs from the time of last scattering to today, so the growth in h g (r) in the bouncing branch at late times 
causes a large increase in <df and therefore in Cf 

to = Hq, then the value of Cj for the bouncing branch is 10 11 times larger than for pure GR: C. 


We find that if we set h g (Ti) = hf(n) at n = 10 6 H 0 1 and 

n T 


2 BB 


10 1 


y 2 GR’ 


demonstrating that the growth at late times dominates the signal. However, our initial time r,; = 10 _b H^ corre¬ 
sponds to a reheat temperature of only Tj = 0.07 GeV. Extrapolating to a more reasonable reheat temperature, say 
Ti = 10 10 GeV, requires looking to equation (291 which shows that h g grows with T,. and so we expect an even bigger 
enhancement of the quadrupole. 


To ensure that Cj BB = C^qj^ ~ 200/iA' 2 as phenomenologically required, we must divide the spectrum by a large 
factor, which is equivalent to a large suppression of the initial overall amplitude of tensor perturbations hig j^Ti). 
This can be accomplished by e.g. lowering the energy scale of inflation. In addition, by tuning the initial conditions, 
we can control the growth in h g , and relieve the need for this large suppression factor. Referring to equation ( |29| , 
we see that an adjustment of hf(Ti)/h g (Ti ), 7$ (or Tj), or m can lead to a power spectrum Cf that is consistent with 
observations. For example, FIG. [5] shows how the bouncing branch power spectrum converges to the standard one as 
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hf(Ti)/h g ( T i) is decreased. Here, we have re-scaled the power spectrum to obtain Cj/Cf ~ 0.2. Note that these plots 
would look the same if instead, m or T) was decreased by the same amount. In summary, to achieve a CMB Tensor 
Power Spectrum that resembles the result from GR, we require very tuned initial conditions, or tuned graviton mass. 





FIG. 5: The CMB Tensor Power Spectrum in the bouncing branch (red, dashed) approaches the GR result (blue, solid) as the initial 
value of hf (or m or T, ) is decreased from left to right: hf (r, ) jh g (r,) = 10 3 , 10 5 , 10 7 with r, = 10 — 6 Hq 1 . In each cases we have 
scaled the power spectrum in the bouncing branch down so that CfT BB ~ 200 [iK 2 , which required dividing the spectrum by 
~ 10 6 , 5.3, 1.1 from left to right. 


V. PRESENT DAY STOCHASTIC GRAVITATIONAL WAVE BACKGROUND 


We now use the results from the previous section to see how the bigravity primordial gravitational waves contribute 
to the present day stochastic gravitational wave energy density . The observable quantity of interest is the gravitational 
wave energy density, defined as a function of frequency: 


^GW (/) 


1 dpow 
Pcrit din f 


(32) 


where the critical density is p cr i t = 3M 2 H 2 (t). 


Direct detection of relic gravitational waves is of considerable interest given the improving technology of ground 
and space based laser interferometers. Various experiments have already placed upper bounds on Bq W , and proposed 
experiments will be able to reach much higher sensitivities. Therefore, one might ask if gravitational waves in bigravity 
would be more or less likely to detect, and if the current sensitives of LIGO or Pulsar Timing Arrays could constrain 
bigravity. LIGO has already made measurements between between 51 < / < 150 Hz to constrain Bq W < 6.5 x 10 -5 
at these frequencies and advanced LIGO is predicted to reach down to sensitivities of Qq W ~ 6.5 x 10~ 9 in the 
coming years [34j . In addition, Pulsar-timing experiments have placed an upper bound of Bq W < 1.6 x 10~ 9 at 
low frequencies 10 -9 < / < 10~ 8 Hz |35j and will improve in the future. The first-generation space based laser 
interferometer, LISA, is expected to operate at sensitivities of 12 q W ~ 10” 11 at frequencies / ~ 10 -3 Hz [3B], while 
the second-generation space based interferometer, BBO, may be able to reach all the way down to Oq W ~ 10 —1 ' near 
frequencies / ~ 0.3 Hz EZI- 


In terms of the tensor modes (corresponding to the physical metric) the predicted stochastic background can be 
computed at any conformal time r via the formula (38j 


^GW (fc,r) 


k 2 \h g (k,T)\ 2 + \h g (k,T)\ 2 

12t t 2 H 2 (t) 


(33) 


given as a function of k = 2irf. When written with superscript 0, it is understood to be evaluated today r = tq. More 
generally however, Hq W (/c) represents the present-day gravitational wave energy density on scales that re-entered the 
Hubble horizon during the radiation dominated era. 


Over the range of frequencies of interest for the experiments above, the stochastic background due to primordial 
tensor modes in GR is essentially flat Hq W (A:) ~ 10 -15 [39]. The precise profile depends on the assumed model of in¬ 
flation that produced the modes, which for us is unimportant as we are just looking for a rough comparison to bigravity. 
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Let us compare njl w (fc) in bigravity and GR. In the bouncing branch, we have observed that the late time growth 
of tensor modes falls of with the square of the frequency, and therefore, using (33), so will SlQ W (fc). This makes 
Hqw(^) harder to detect in the bouncing branch compared to GR over the range of frequencies of interest for the 
experiments listed above, 10~ 9 < / < 10 3 Hz (or 10 9 < k/H a < 10 21 ). In the expanding branch, the decay of the 
tensor modes closely matches with GR, so in this case we expect a result similar to the standard picture. See FIG. [6] 
for a plot of the results for Hq W over a range of frequencies from k = IOHq to k = 10 4 Ho- 


Bouncing Branch 



Expanding Branch 



FIG. 6: The present day stochastic gravitational wave background, given by \33\ for bigravity as compared to GR. The expanding 
branch (green) agrees closely with GR (blue), while the bouncing branch (red) shows drastic differences. Note that we have scaled £2 q W 
down by an appropriate factor so as to fix C^ BB ~ 200^i^ 2 (see FIG. pl. 


As in the previous section, we see that an adjustment of the initial condition for hf causes the result for the 
bouncing branch to converge to the solution in GR, as displayed in FIG. [7] This is equivalent to varying m or Ti 
by the same amount, as discussed previously. In this plot, we have re-scaled the power spectrum to yield a tensor 
contribution to the CMB temperature quadrupole of Cj BB ~ 200 fiK 2 . 


§ 


CD 


a 



FIG. 7: The present day stochastic gravitational wave background in the bouncing branch (red) approaches the GR result (blue) as the 
initial value of hf (or m or T,) is decreased from left to right: hj (r,)/ hg (r ? ) = 10 3 , 10 —5 , 10 —7 with n = I0 —6 Hq 1 


VI. INITIAL CONDITIONS 

We have seen that in the bouncing branch, the extreme growth in the dark sector causes amplification of the 
physical tensor mode, leading to large discrepancies with GR. This amplification causes alterations in physical ob¬ 
servables, such as the CMB Power Spectrum and the present day stochastic gravitational wave background. However, 
if the tensor modes in the dark sector are sufficiently suppressed, then the physical tensor modes and their associ¬ 
ated observables closely resemble those of GR. If some mechanism were to exist so that /i/(t,;) <C h g (ri) then this 
branch would have essentially the same gravitational wave spectrum as GR, and would be indistinguishable on an 
observational level. The need for this tuning has been observed in nn nz], and further in EZ|, where they estimate 












12 


the required suppression to match the observed CMB power spectrum. It is therefore necessary to explore the initial 
conditions for the primordial tensor modes, assuming they were produced by inflation. 


During inflation, the universe undergoes accelerated expansion in a quasi-de Sitter phase. For pure de Sitter, the 
bigravity background equations simplify as follows (see also Ref. [40]): 


p = constant = 3(HgM g ) 2 => p = 3(Hg/H 0 ) 2 
r = constant 

c= 1 => n = H g = 'Hf 


(34) 

(35) 

(36) 


where H 1 is the Hubble parameter during inflation. 


The second line follows from equation (22) and the third line 


follows from (13). Therefore, the dark universe is also undergoing de Sitter expansion. Restoring m and Hq in (22) 


using (14), specializing to the bouncing branch parameters in which only ft and ft are nonzero, we get a polynomial 
equation for the value of the ratio of the scale factors during inflation, r 1 : 


- 3(ff 7 )V + m 2 [ft + ft(r 7 ) 3 - 3ft(ft) 2 ] = 0. 


(37) 


Since one typically takes the mass term rn in bigravity to be on the order of the Hubble constant ift which is much 
smaller than the Hubble parameter during inflation H g , we can expand the solution for r 1 in powers of Hg/m 3 > 1 , 
yielding 


r 


/ 



Hi 

Hn 


3ft 
2ft 




(38) 

(39) 


where we have used m = Hq => ft = ft = 0.94. Assuming high scale inflation, the maximum allowed Hubble during 
inflation is Hg ~ 10 15 GeV, and Hq ~ 10 -33 eV, from which we can estimate ft ~ 10 57 . For very low scale inflation, 
say at the TeV scale, we can estimate ft ~ 10 29 . 


After inflation, there must be a transition to a radiation dominated phase. During the inflationary de Sitter phase, 
the scale factor for the dark sector metric is increasing, b > 0, but in the radiation dominated era in the bouncing 
branch, the scale factor is decreasing, b < 0. It is evident that the dark sector must undergo another “bounce” tran¬ 
sition after inflation from expansion to contraction in order to achieve the necessary behaviour in the early radiation 
era. In the next section we see that this bounce is indeed achieved in a simple inflationary model. As argued for the 
late-tinre bounce in this branch, it is likely that the associated divergence in the curvature is a mathematical feature 
rather than a physical problem. 


A note on inflation in the expanding branch is now in order. In this branch, we set ft 
with only ft ^ 0, we obtain 


0. Solving for r 1 in (37) 


r 


i 


-3 Hf + y/9H* 2 + 12 m 4 ft 2 
6?T7 2 ft 


ftm 2 

3ftP 


Ho 

Hl g 


O 



Therefore, in the expanding branch we obtain a very small value of r 1 . 


(40) 

(41) 

(42) 


To find the power spectrum of primordial tensors, we expand the action Eq. 10 to quadratic order in transverse 
traceless perturbations of the / and g metrics. During inflation, the interaction terms are unimportant due to the 
large hierarchy between m and Hg. Defining the canonically normalized fields: 


Vn = 


MgCL 

—3—h n 


v f = 




2 


2 


(43) 













13 


and using the fact that 

a b 2 
a b t 2 

in de Sitter, the quadratic action for tensors during an inflationary epoch in bigravity is given by: 


5 = E \ I dr d3k 


+,X 




u gJ 


Imposing Bunch Davies initial conditions, the mode functions each obey: 

1 


%/ = 


k 3 / 2 : 


(1 — ikr) e 


ikr 


Summing over the two polarization states, the power spectrum for hf and h g are given by: 


= 


2k 3 


7 t 2 M 2 a 2 


,(/) 2 _ 2fc 3 \v f \' 2 


k=H 


A U ’ = 

T n 2 M 2 b 2 


k=U 


Substituting with Eqn. 46 and using H g = (ar) 1 for a de Sitter phase, we obtain: 


( 44 ) 


(45) 


(46) 


(47) 


A (s ) 2 = 2H L 

T 7 T 2 M 2 


2 

a (/) 2 _ A T 

— t 2 


k=H 


(48) 


Modes are populated for both the visible and dark sector tensors on all super horizon scales. The appearance of 
r 1 is a consequence of the fact that the relative size of a and b is physical, and cannot be removed by a change of 
coordinates or field redefinition. 


In the bouncing branch r 1 1, leading to a drastic suppression in the initial amplitude for hf roughly given 
by hf/hg = 1 /r 1 ~ Hq/Hq. This suppression is more than sufficient to bring the amplitude of dark sector tensor 
fluctuations below the threshold where they alter the propagation of visible sector tensors. With inflationary initial 
conditions, we therefore conclude that there would be no visible deviation from GR in the tensor contribution to the 
CMB temperature anisotropies or the late time stochastic distribution of gravitational waves. This is true for both 
high scale or low scale inflation. 

In contrast, the expanding branch gives hf/hg ~ (Hg /Hq) 2 , which for any reasonable choice of the inflationary 
scale is far beyond the perturbative regime. Therefore, we cannot make sense of inflation in the expanding branch 
with /3 4 = 0. If instead we considered a non-minimal expanding branch with 04 ^ 0, then the result would be the 
same as the bouncing branch. 


Finally, let us comment on the time dependence of the dark sector tensor modes. This is of interest because as 
illustrated in mm, the dark sector tensor mode has a growing mode on super horizon scales proportional to r 3 
during radiation domination. We can estimate the relative amplitude of the growing and constant modes of (28) at 
the beginning of radiation domination as |r/i/(Tj)//i/(Tj)|. Using (43) and (46), one obtains 


IMu)l = 


2 (1 + k 2 r 2 ) 1 / 2 2 H^l + k 2 ^) 1 / 2 


M g k 3 / 2 Tib(Ti) M g 


fc 3 / 2 r J 


IM T i)| = 


Hgk 2 Ti 


M g 2k 3 / 2 r I (l + k 2 T 2 y/ 2 ’ 


(49) 


where we used that 6(r*) = r / a(r;) = r 1 / (HgTi). In the limit of small iy, the ratio between the growing and constant 
mode is therefore 


Tihf(n) 


hf( T i) 


2 a(Ti)m- 


Using a ex t during radiation domination, for horizon-scale wave numbers k 
estimate this ratio as 


Tihf(Ti) 


hf(n) 


HnTi 


i/ T o)Hl 


^0 


10 


-57 


(50) 

Hq, in the bouncing branch we can 

(51) 
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For any reasonable choice of reheat temperature, this will be smaller than the growth factor for the growing mode 
during radiation domination, ~ (T eq /ri ) 3 = (Treh/Teq ) 3 ~ 10 30 , where the subscript “eq” refers to matter radiation 
equality. We conclude that the growing mode remains significantly suppressed, and therefore the appropriate initial 
conditions are h g j = 0 . 

Before moving on to a specific inflationary model in bigravity, let us comment on the validity of these calculations. 
Massive bigravity is an effective field theory with a strong coupling scale A 3 = ( m 2 M g ) 1 / 3 . For m ~ H 0 , this scale 
is rather low: A 3 ~ 10 -22 GeV. This is much smaller than any reasonable inflationary scale. Therefore, we should 
worry about the validity of the effective field theory during inflation. While the scalar perturbations become strongly 
coupled, the potential term for the transverse traceless modes becomes irrelevant during inflation. However, we should 
be concerned about interactions between the scalar and transverse traceless perturbations which may contribute to 
the evolution of the tensors. A possible resolution to describing inflation in bigravity lies in the Vainshtein screening 
mechanism m- If screening is efficient in the high density environment of inflation, the scale strong coupling scale 
gets redressed to a far higher scale A| A 3 , and the structure of the potential is not affected (see, for instance 
®;43j). Under this assumption of screening, calculations at the inflationary scale might still be valid in the effective 
field theory description of bigravity. A rigorous proof of this argument is still required, which we leave to future work. 


A. An inflationary model in bigravity 


In this section, we examine bigravity for the m 2 (f > 2 inflationary model. We want to determine the behaviour of 
r (and thus b) during inflation in the bouncing branch with only 4 ^ 0. The inflaton field has potential energy 
V(<f>) = ^m^cj) 2 and energy density p$ = \(d t (j )) 2 + V(<f>). For this calculation we find it convenient to use the following 
set of dimensionless variables: 


H 

H -V 




t = tHi 


T rb = V = 


V 




R = — _ R* = g 0 R 

Pn ( H iyPn 


Pr — 


Mg (Hi ) 2 
Pr 

M 2 (my 


(52) 

(53) 


where is the decay rate of the inflaton, t is proper time, and we take (H g ) 2 = V((f>o) /3M 2 in terms of the value of 
</> at the start of inflation, implying that V = 3 cp 2 /<%. 


Including an explicit decay of the inflaton into radiation, the equation of motion is 


4>" + 3 H# + + 6 A = 0 

<t> o 


(54) 


where a prime denotes a derivative with respect to dimensionless proper time t. 
a(t ) becomes 


The Friedmann equation (151 for 


3 H 2 — Pfy + p r — ~((j )') 2 + 3 i- + p r - (55) 

2 n 

Notice that we have neglected the contribution of p m and p mg in the early universe since these will be highly suppressed 
compared to the inflaton or radiation energy density. The radiation energy density p r satisfies a modified conservation 
equation 

Pr + AHp r = f (56) 


We can now solve (54), ( |55| ), and ( 56J) for the functions </>(t), a(i), p r (t). The last ingredient will be to solve for r(f), 
for which we use ( 22 ), which simpl: 


ifies in the bouncing branch with /3-| .4 ^0 to 
0 = rp — fix — jj^r 3 + 3/3ir 2 

H l IT r , 

r = -f-\l — tor large r 

tin \ Da 


(57) 

(58) 


The results of the calculation are shown in FIG. [ 8 ] We can see the transition from b 1 > 0 during inflation to b 1 < 0 
after inflation is achieved in the bottom right plot. Notice that after b hits its first peak, it oscillates as it decreases, 
indicating brief periods of expansion and contraction of This behaviour is caused by the oscillation of the inflaton 
around its minimum, and implies that the dark sector metric undergoes multiple bounces during reheating. 
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FIG. 8: Top left: The solutions of ( |54[ ), ( |55[ |, and ( |56| ) for </>, p <£, and p r vs. t for the m 2 (j) 2 model of inflation. We set T$ = 10 3 , 
0(0) = </>o = 24, 0'(O) = 0, and p r ( 0) = 0. Top right: The scale factor for Bottom left: The ratio r found via ( |58| >. Bottom 
right: The scale factor for b = ra. 


B. The effects of reheating on tensor perturbations 


To determine the evolution of perturbations during the transition from inflation to a radiation dominated Universe, 
we have explicitly evolved the coupled evolution equations Eq. (251 and (261 in the background shown in Fig. [8j The 
evolution of a few modes from inflation, through reheating, to near the end of radiation domination, with initial con¬ 
ditions set by Eq. (461, is shown in Fig. [9j On scales relevant for the CMB, we find that after horizon-crossing during 
inflation hf remains frozen through reheating and the duration of radiation domination. This validates the initial 
conditions assumed above, where only the constant mode of hf is excited at the beginning of radiation domination. 
On these large scales, we find good quantitative agreement with the prediction Eq. (51), and conclude inflationary 
initial conditions will yield no observable deviation from GR in the CMB. 


On far smaller scales, of order the size of the horizon at reheating (in our model, this corresponds to k ~ 10 3O Ho), 
the growing mode is excited early in the radiation dominated era. For such modes, one should not apply the growth 
factor for h g found above (Eq. ([29])) , but rather the growth factor found in Ref. H7]. Since this growth factor falls off 
like fc 4 , the growth in h g is highly suppressed, and the result will be no observable deviation from GR even though 
the growing mode dominates from the beginning of radiation domination. 


VII. CONCLUSION 

We have studied the properties of gravitational waves in massive bigravity, and their impact on cosmological 
observables compared to the standard predictions of General Relativity. The two background solutions we have stud¬ 
ied display dramatically different phenomenology, illustrating the enormous size of the parameter space for observables. 

In the “expanding branch” in which both metrics expand in time, we found that for a wide range of initial condi- 
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FIG. 9: The evolution of hf from inflation (at redshift of z ~ 10 54 ), through reheating (at redshift z ~ 10 30 ), to the end of radiation 
domination (at redshift z ~ 10 3 ). 


tions, the physical tensor perturbations h g matches closely with the pure GR solution. Due to a dramatic decay of 
hf , the impact of the dark sector on h g is not important and causes no significant deviation from GR. 

The “bouncing branch”, in which the dark metric f gl , undergoes a bounce from contraction to expansion, has 
potentially dramatic differences from GR in the tensor sector. When the / metric is undergoing contraction, the lapse 
c is negative, which causes hf to grow. This growth in hf can translate into growth in h g through the mixing term in 
the equations of motion, in some cases leading to physical gravitational waves with growing amplitudes at late times. 
This contrasts significantly with gravitational waves in GR which decay with time. The growth can potentially impact 
the CMB tensor power spectrum by dramatically amplifying large scale temperature anisotropies. The present day 
stochastic gravitational wave background, Hq W can be impacted as well, inheriting a very red spectrum that decays 
with the square of the wave number. The degree of growth depends on the scale of reheating, amplitude of the dark 
sector tensor modes, wave number, and graviton mass, and obeys the scaling relation Eq. [29] for initial amplitudes 
above a critical value. On the largest scales, we find that the dark sector tensor modes can have a significant influence 
on the physical tensor for hf(Ti)/h g (Ti) > 10 -20 for a reheat temperature of T* = 10 10 GeV. In the absence of a 
theory of initial conditions, it is not clear that this holds. 

To address the question of initial conditions, we computed the primordial power spectrum for dark and visible sector 
tensors in an inflationary cosmology. We found that the expanding branch is far beyond the perturbative regime, and 
therefore inaccessible to a semi-classical treatment. However, the primordial power spectra in the bouncing branch 
show that hf/h g ~ H 0 /Hg ~ 10 -57 for high scale inflation, and hf/h g ~ 10~ 29 for low scale inflation. With this 
level of suppression, there will be no observable deviation from GR in the CMB or stochastic gravitational wave 
background. We presented an inflationary model that exhibits this explicitly. 

Let us now discuss our results in the context of related work in the literature. While this work was in progress, 
Refs, mm appeared with complementary investigations of tensor modes in bigravity. Ref. m considered the 
bouncing branch chosen in this paper with identical parameters, but with an initial condition that was entirely 
composed of a growing mode of the dark sector tensor. However, as shown in Sec. |VI| it appears that inflationary 
initial conditions do not excite the growing mode as significantly. In m the authors considered varying the initial 
conditions in a phenomenological way, and specifically tested the effects on CMB spectra. Both of these works 
conclude that a tuning of the initial conditions can render the theory viable, which we show is actually possible in the 
inflationary paradigm. In this sense, our investigation is largely complementary to previous work; taken together, the 
range of possibilities is covered. Of course, another theory of initial conditions may prevail, and a proper treatment 
of higher-order couplings between the scalar and tensor sectors may reveal a significant enhancement. However, in 
the context of linear perturbation theory in inflationary cosmology, it appears that the growing mode on large scales 
is not excited at the end of inflation. 


It is clear that tensors can be a sensitive probe of massive bigravity. Looking to the future, the parameter space 
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of nearly homogeneous solutions will soon be completely explored both at the level of the background and first order 
perturbations. In light of this, it is equally important to consider the theory for initial conditions in a broader sense, 
as illustrated by the strong dependence on initial conditions found in this and other papers. To this end, we plan 
to return to the question of inflationary model building in massive bigravity in future work. Scenarios with small 
but observable deviations from GR could serve as an important alternative hypothesis necessary for testing GR on 
cosmological scales and in future gravitational wave observatories. 
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